#include <stdlib.h>

#include "../image.h"
#include <util/math_macros.h>
#include <util/msg.h>

/* ========================================================================== */

#define LABEL_INC 10;

/* ========================================================================== */
// <editor-fold defaultstate="collapsed" desc="Auxiliary Functions to label_he_chao_suzuki_wu">

inline int * __create_table(int size) {
	int i;
	int * table = mem_allocn(size, int);
	for (i = 0; i < size; i++) table[i] = -1;
	return table;
}

inline static int __label_find(int x, int * Parent) {
	while (Parent[x] != -1) x = Parent[x];
	return x;

}

inline static void __label_union(int x1, int x2, int * Parent) {
	x1 = __label_find(x1, Parent);
	x2 = __label_find(x2, Parent);
	if (x1 != x2) Parent[x2] = x1;
}
// </editor-fold>

/***************************************************************************
 * This function is implemented based on work:
 * 
 *	Lifeng He, Yuyan Chao, Kenji Suzuki, and Kesheng Wu. 2009. 
 *	Fast connected-component labeling. Pattern Recogn. 42, 9
 *	(September 2009), 1977-1987. DOI=10.1016/j.patcog.2008.10.013
 *	http://dx.doi.org/10.1016/j.patcog.2008.10.013
 * 
 * @param rows
 * @param columns
 * @param mat
 * @param label
 * @param count
 * @param bv
 * @return 
 */
int label_he_chao_suzuki_wu(int rows, int columns, real * mat, real * label, real fv) {

	int msize;
	int max_comp;
	int nclass;

	int r, c, _r, i;
	int * Table;

	msize = rows*columns;

	max_comp = CEIL(rows / 2.0) * CEIL(columns / 2.0);

	Table = __create_table(max_comp);

	{/* Step 01: Forward raster scan ======================================= */
		real left;
		real top_left;
		real top;
		real top_right;
		real c1, c2, c3, c4;
		real LABEL = 1;

		nclass = 0;

		for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
			for (c = 1; c < columns - 1; c++) {
				i = _r + c;

				if (mat[i] != fv) continue;

				/* = mat[i] is an object ==================================== */

				left = label[i - 1];
				top = label[i - columns];
				top_left = label[ i - columns - 1];
				top_right = label[i - columns + 1];

				c1 = left;
				c2 = top_left;
				c3 = top;
				c4 = top_right;

				if (c3 != 0) label[i] = c3;
				else if (c1 != 0) {
					label[i] = c1;
					if (c4 != 0) __label_union(c4, c1, Table);
				} else if (c2 != 0) {
					label[i] = c2;
					if (c4 != 0) __label_union(c2, c4, Table);
				} else if (c4) label[i] = c4;
				else {
					// new label
					label[i] = LABEL;
					LABEL += 1;
					nclass++;
				}
			}
		}
	}
	{ /* Step 01.5: Intermediate step to compress the parent table ========== */
		for (i = 0; i < max_comp; i++) {
			if (Table[i] == -1) continue;
			nclass--;
			Table[i] = __label_find(i, Table);
		}
	}
	{ /* Step 02: Backward verification ===================================== */
		long int x;
		for (i = 0; i < msize; i++) {
			x = label[i];
			label[i] = __label_find(x, Table);
		}
	}
	{ /* Step 03: Convert all labels to sequence 0, 1, 2, ... , *count ====== */
		int c;
		int p;
		c = 1;
		for (i = 0; i < msize; i++) {

			p = (int) label[ i ];

			/* Is background */
			if (p == 0) continue;

			/* Is foreground */
			if (Table[ p ] == -1) {
				Table[ p ] = c;
				c++;
			}
			label[ i ] = Table[ p ];
		}
	}
	free(Table);

	return nclass;
}

int imp_image_label(ImpMatrix * Mat, ImpMatrix ** Label, real fv) {

	int nclass;
	if (*Label == NULL) (*Label) = imp_matrix_new(Mat->rows, Mat->columns);

	nclass = label_he_chao_suzuki_wu(Mat->rows, Mat->columns, Mat->data, (*Label)->data, fv);

	return nclass;
}

void __imp_image_label(ImpMatrix * Mat, ImpMatrix ** Label, int * count, real VF) {

	real * mat;
	real * label;

	int rows, columns;
	long int msize;

	long int max_comp;

	int * Table;

	int r, c;
	long int _r, i;

	_verify(Mat);
	_verify(Label);

	rows = Mat->rows;
	columns = Mat->columns;
	msize = Mat->size;
	mat = Mat->data;

	if (*Label == NULL) (*Label) = imp_matrix_new(rows, columns);
	label = (*Label)->data;

	max_comp = CEIL(rows / 2.0) * CEIL(columns / 2.0);

	Table = (int *) calloc(max_comp, sizeof (int));
	for (i = 0; i < max_comp; i++) {
		Table[i] = -1;
	}

	{ /* Step 01: Forward raster scan ======================================= */
		real left;
		real top_left;
		real top;
		real top_right;

		real c1, c2, c3, c4;

		real LABEL = 1;

		*count = 0;

		for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
			for (c = 1; c < columns - 1; c++) {
				i = _r + c;

				if (mat[i] != VF) continue;

				/* = mat[i] is an object ==================================== */

				left = label[i - 1];
				top = label[i - columns];
				top_left = label[ i - columns - 1];
				top_right = label[i - columns + 1];

				c1 = left;
				c2 = top_left;
				c3 = top;
				c4 = top_right;

				if (c3 != 0) label[i] = c3;
				else if (c1 != 0) {
					label[i] = c1;
					if (c4 != 0) __label_union(c4, c1, Table);
				} else if (c2 != 0) {
					label[i] = c2;
					if (c4 != 0) __label_union(c2, c4, Table);
				} else if (c4) label[i] = c4;
				else {
					// new label
					label[i] = LABEL;
					LABEL += 1;
					(*count)++;
				}
			}
		}
	}
	{ /* Step 01.5: Intermediate step to compress the parent table ========== */

		for (i = 0; i < max_comp; i++) {
			if (Table[i] == -1) continue;

			(*count)--;
			Table[i] = __label_find(i, Table);
		}

	}
	{ /* Step 02: Backward verification ===================================== */
		long int x;
		for (i = 0; i < msize; i++) {
			x = label[i];
			label[i] = __label_find(x, Table);
		}
	}
	{ /* Step 03: Convert all labels to sequence 0, 1, 2, ... , *count ====== */
		long int c;
		long int p;
		c = 1;
		for (i = 0; i < msize; i++) {

			p = (long int) label[ i ];

			/* Is background */
			if (p == 0) continue;

			/* Is foreground */
			if (Table[ p ] == -1) {
				Table[ p ] = c;
				c++;
			}
			label[ i ] = Table[ p ];
		}
	}

	free(Table);
}

void imp_image_label_analisys(ImpMatrix * Label, int nclass, ImpMatrix ** Feat) {

	enum _features {
		X0 = 0,
		Y0,
		X1,
		Y1,
		Xc,
		Yc,
		AREA,
		RATIO,
		Xcc,
		Ycc,
		N_FEATURES
	};

	int r, c, _r, i, f;
	int cl;
	int rows, columns;
	real * feat;
	real * label;

	*Feat = imp_matrix_new(nclass, N_FEATURES);
	feat = (*Feat)->data;

	label = Label->data;

	rows = Label->rows;
	columns = Label->columns;

	// = Initialization
	for (r = 0; r < nclass; r++) {

		f = r * N_FEATURES;
		feat[f + X0] = Label->size;
		feat[f + Y0] = Label->size;
		feat[f + X1] = 0;
		feat[f + Y1] = 0;
		feat[f + Xc] = 0;
		feat[f + Yc] = 0;
		feat[f + AREA] = 0;
		feat[f + RATIO] = 0;
		feat[f + Xcc] = 0;
		feat[f + Ycc] = 0;
	}

	for (r = 0, _r = r * columns; r < rows; r++, _r += columns) {
		for (c = 0; c < Label->columns; c++) {
			i = _r + c;

			cl = (int) label[i];

			if (cl == 0) continue;

			cl = cl - 1;

			f = cl * N_FEATURES;

			if (feat[f + X0] > c) feat[f + X0] = c;
			if (feat[f + Y0] > r) feat[f + Y0] = r;

			if (feat[f + X1] < c) feat[f + X1] = c;
			if (feat[f + Y1] < r) feat[f + Y1] = r;

			feat[f + AREA]++;

			feat[ f + Xcc ] += c;
			feat[ f + Ycc ] += r;

		}
	}

	for (r = 0; r < nclass; r++) {
		f = r * N_FEATURES;
		feat[ f + Xc ] = (feat[f + X0] + feat[f + X1]) / 2;
		feat[ f + Yc ] = (feat[f + Y0] + feat[f + Y1]) / 2;

		feat[ f + RATIO] = (feat[ f + X1] - feat[ f + X0]) / (feat[ f + Y1] - feat[ f + Y0]);

		feat[ f + Xcc ] = feat[f + Xcc] / (feat[f + AREA]);
		feat[ f + Ycc ] = feat[f + Ycc] / (feat[f + AREA]);
	}



}

/* ========================================================================== */

void old__imp_image_label(ImpMatrix * Mat, ImpMatrix ** Label, int * count) {

	/***************************************************************************
	 * This procedure is implemented based on work:
	 * 
	 *		Lifeng He, Yuyan Chao, Kenji Suzuki, and Kesheng Wu. 2009. 
	 *		Fast connected-component labeling. Pattern Recogn. 42, 9
	 *		(September 2009), 1977-1987. DOI=10.1016/j.patcog.2008.10.013
	 *		http://dx.doi.org/10.1016/j.patcog.2008.10.013
	 * 
	 **************************************************************************/

	real * mat;
	real * label;

	int rows, columns;
	long int msize;

	long int max_comp;

	int * Table;

	int r, c;
	long int _r, i;

	_verify(Mat);
	_verify(Label);

	rows = Mat->rows;
	columns = Mat->columns;
	msize = Mat->size;
	mat = Mat->data;

	if (*Label == NULL) (*Label) = imp_matrix_new(rows, columns);
	label = (*Label)->data;

	max_comp = CEIL(rows / 2.0) * CEIL(columns / 2.0);

	Table = (int *) calloc(max_comp, sizeof (int));
	for (i = 0; i < max_comp; i++) {
		Table[i] = -1;
	}

	{ /* Step 01: Forward raster scan ======================================= */
		real left;
		real top_left;
		real top;
		real top_right;

		real c1, c2, c3, c4;

		real LABEL = 1;

		real VB;
		//real VF;

		VB = 0;
		//VF = 255;

		*count = 0;

		for (r = 1, _r = r * columns; r < rows - 1; r++, _r += columns) {
			for (c = 1; c < columns - 1; c++) {
				i = _r + c;

				if (mat[i] == VB) continue;

				left = label[i - 1];
				top = label[i - columns];
				top_left = label[ i - columns - 1];
				top_right = label[i - columns + 1];

				c1 = left;
				c2 = top_left;
				c3 = top;
				c4 = top_right;

				if (c3 != VB) label[i] = c3;
				else if (c1 != VB) {
					label[i] = c1;
					if (c4 != VB) __label_union(c4, c1, Table);
				} else if (c2 != VB) {
					label[i] = c2;
					if (c4 != VB) __label_union(c2, c4, Table);
				} else if (c4) label[i] = c4;
				else {
					// new label
					label[i] = LABEL;
					LABEL += 1;
					(*count)++;
				}
			}
		}
	}
	{ /* Step 01.5: Intermediate step to compress the parent table ========== */

		for (i = 0; i < max_comp; i++) {
			if (Table[i] == -1) continue;

			(*count)--;
			Table[i] = __label_find(i, Table);
		}

	}
	{ /* Step 02: Backward verification ===================================== */
		long int x;
		for (i = 0; i < msize; i++) {
			x = label[i];
			label[i] = __label_find(x, Table);
		}
	}

	free(Table);
}

void _select8(ImpImage Gray, int x, int y, ImpImage ** Sel, unsigned char px) {

	unsigned char * gray;
	unsigned char * sel;

	int rows;
	int columns;

	int i;

	rows = Gray.height;
	columns = Gray.width;

	if (x < 0 || x >= columns) return;
	if (y < 0 || y >= rows) return;

	_verify(Sel);
	if (Sel) if (*Sel == NULL) (*Sel) = imp_image_new(Gray.width, Gray.height, IMP_IMAGE_TYPE_GRAY);


	gray = Gray.data;
	sel = (*Sel)->data;

	i = y * columns + x;

	if (sel[i] > 0) {
		// pixel just visited
		return;
	}

	if (gray[i] == 0) {
		// not connected pixel
		//sel[i] = 100; // Mark border
		return;
	}



	sel[i] = px;

	// left pixel
	_select8(Gray, x - 1, y, Sel, px);

	// right pixel
	_select8(Gray, x + 1, y, Sel, px);

	// top pixel
	_select8(Gray, x, y - 1, Sel, px);

	// bottom pixel
	_select8(Gray, x, y + 1, Sel, px);

	// top left
	_select8(Gray, x - 1, y - 1, Sel, px);

	// top right
	_select8(Gray, x + 1, y - 1, Sel, px);
	// bottom left
	_select8(Gray, x - 1, y + 1, Sel, px);
	// bottom right
	_select8(Gray, x + 1, y + 1, Sel, px);


}

void _select4(ImpImage Gray, int x, int y, ImpImage ** Sel, unsigned char px) {

	unsigned char * gray;
	unsigned char * sel;

	int rows;
	int columns;

	int i;

	rows = Gray.height;
	columns = Gray.width;

	if (x < 0 || x >= columns) return;
	if (y < 0 || y >= rows) return;

	_verify(Sel);
	if (Sel) if (*Sel == NULL) (*Sel) = imp_image_new(Gray.width, Gray.height, IMP_IMAGE_TYPE_GRAY);


	gray = Gray.data;
	sel = (*Sel)->data;

	i = y * columns + x;

	if (sel[i] > 0) {
		// pixel just visited
		return;
	}

	if (gray[i] == 0) {
		// not connected pixel
		//sel[i] = 100; // Mark border
		return;
	}



	sel[i] = px;

	// left pixel
	_select4(Gray, x - 1, y, Sel, px);

	// right pixel
	_select4(Gray, x + 1, y, Sel, px);

	// top pixel
	_select4(Gray, x, y - 1, Sel, px);

	// bottom pixel
	_select4(Gray, x, y + 1, Sel, px);

}

/*
void _old_imp_image_label8(ImpImage Gray, ImpImage ** Label, int * count) {

	unsigned char * gray;
	unsigned char * label;

	int rows;
	int columns;

	int r, c, _r, i;
	int up, left;

	unsigned char LABEL;

	_verify(Label);
	if (Label) if (*Label == NULL) (*Label) = imp_image_new(Gray.width, Gray.height, IMP_IMAGE_TYPE_GRAY);

	gray = Gray.data;
	rows = Gray.height;
	columns = Gray.width;

	label = (*Label)->data;

	LABEL = 0;
	(*count) = 0;

	// Step 01:
	for (r = 0, _r = r * columns; r < rows; r++, _r += columns) {
		for (c = 0; c < columns; c++) {
			i = _r + c;

			if (gray[i] == 0) continue;

			if (label[i] > 0) continue;

			up = i - columns;
			left = i - 1;

			LABEL += 50;
			(*count)++;

			_select8(Gray, c, r, Label, LABEL);

		}
	}
}

void _old_imp_image_label4(ImpImage Gray, ImpImage ** Label, int * count) {

	unsigned char * gray;
	unsigned char * label;

	int rows;
	int columns;

	int r, c, _r, i;
	int up, left;

	unsigned char LABEL;

	_verify(Label);
	if (Label) if (*Label == NULL) (*Label) = imp_image_new(Gray.width, Gray.height, IMP_IMAGE_TYPE_GRAY);

	gray = Gray.data;
	rows = Gray.height;
	columns = Gray.width;

	label = (*Label)->data;

	LABEL = 0;
	(*count) = 0;

	// Step 01:
	for (r = 0, _r = r * columns; r < rows; r++, _r += columns) {
		for (c = 0; c < columns; c++) {
			i = _r + c;

			if (gray[i] == 0) continue;

			if (label[i] > 0) continue;

			up = i - columns;
			left = i - 1;

			LABEL += 50;
			(*count)++;

			_select4(Gray, c, r, Label, LABEL);

		}
	}
}
 */

